subroutine rks3
use var
implicit none

    integer i,m,j
    real(8) du,dv
    real(8) Qmean_n_k(1:gridsx,4)
    
    
    !k=1
    
    call comput_flux1
    !call comput_flux2    
    call rhs(Qmean,Qmean_n_k)
    do i=1,gridsx
        !write(*,'(I4,2X,100(ES12.4,1X))')i,flux1(i,2)
        Qmean(i,1)=Qmean_n(i,1)-dt*Qmean_n_k(i,1)
            
        Qmean(i,2)=Qmean_n(i,2)-dt*Qmean_n_k(i,2)
            
        Qmean(i,3)=Qmean_n(i,3)-dt*Qmean_n_k(i,3)
            
        Qmean(i,4)=Qmean_n(i,4)-dt*Qmean_n_k(i,4)
        !Qmean(i,j,4)=Qmean_n(i,j,4)-dt*Qmean_n_k(i,j,4)
            
        !write(*,'(I4,2X,100(ES12.4,1X))')i,j,Qmean(i,j,4)
        !do m=1,4
        !    du=-(flux1(i,j,m)-flux1(i-1,j,m))/dx
        !    dv=-(flux2(i,j,m)-flux2(i,j-1,m))/dy
        !    Qmean(i,j,m)=Qmean_n(i,j,m)+dt*du+dt*dv
        !end do            
        !Qmean(i,j,5)=Qmean_n(i,j,5)-dt*Qmean_n_k(i,j,5)        
    enddo 
    call comput_original_var
    call boundary
    
    !k=2
    call comput_flux1
    !call comput_flux2    
    call rhs(Qmean,Qmean_n_k)
     do i=1,gridsx
        Qmean(i,1)=3.d0/4.d0*Qmean_n(i,1)+1.d0/4.d0*(Qmean(i,1)-dt*Qmean_n_k(i,1))!(flux1(i,1)-flux1(i-1,1))/dx
            
        Qmean(i,2)=3.d0/4.d0*Qmean_n(i,2)+1.d0/4.d0*(Qmean(i,2)-dt*Qmean_n_k(i,2))!(flux1(i,2)-flux1(i-1,2))/dx
    
        Qmean(i,3)=3.d0/4.d0*Qmean_n(i,3)+1.d0/4.d0*(Qmean(i,3)-dt*Qmean_n_k(i,3))!(flux1(i,3)-flux1(i-1,3))/dx
            
        Qmean(i,4)=3.d0/4.d0*Qmean_n(i,4)+1.d0/4.d0*(Qmean(i,4)-dt*Qmean_n_k(i,4))!(flux1(i,3)-flux1(i-1,3))/dx
        !Qmean(i,j,4)=3.d0/4.d0*Qmean_n(i,j,4)+1.d0/4.d0*(Qmean(i,j,4)-dt*Qmean_n_k(i,j,4))
        !do m=1,4
        !    du=-(flux1(i,j,m)-flux1(i-1,j,m))/dx
        !    dv=-(flux2(i,j,m)-flux2(i,j-1,m))/dy
        !    Qmean(i,j,m)=3.d0/4.d0*Qmean_n(i,j,m)+1.d0/4.d0*(Qmean(i,j,m)+dt*du+dt*dv)
        !end do
        !Qmean(i,j,5)=3.d0/4.d0*Qmean_n(i,j,5)+1.d0/4.d0*(Qmean(i,j,5)-dt*Qmean_n_k(i,j,5))       
    enddo
    call comput_original_var
    call boundary
    
    !k=3
    call comput_flux1
    !call comput_flux2
    call rhs(Qmean,Qmean_n_k)
    do i=1,gridsx
        Qmean(i,1)=1.d0/3.d0*Qmean_n(i,1)+2.d0/3.d0*(Qmean(i,1)-dt*Qmean_n_k(i,1))!(flux1(i,1)-flux1(i-1,1))/dx
        
        Qmean(i,2)=1.d0/3.d0*Qmean_n(i,2)+2.d0/3.d0*(Qmean(i,2)-dt*Qmean_n_k(i,2))!(flux1(i,2)-flux1(i-1,2))/dx
            
        Qmean(i,3)=1.d0/3.d0*Qmean_n(i,3)+2.d0/3.d0*(Qmean(i,3)-dt*Qmean_n_k(i,3))!(flux1(i,3)-flux1(i-1,3))/dx
            
        Qmean(i,4)=1.d0/3.d0*Qmean_n(i,4)+2.d0/3.d0*(Qmean(i,4)-dt*Qmean_n_k(i,4))
        !do m=1,4
        !    du=-(flux1(i,j,m)-flux1(i-1,j,m))/dx
        !    dv=-(flux2(i,j,m)-flux2(i,j-1,m))/dy
        !    Qmean(i,j,m)=1.d0/3.d0*Qmean_n(i,j,m)+2.d0/3.d0*(Qmean(i,j,m)+dt*du+dt*dv)
        !end do
        !Qmean(i,j,5)=1.d0/3.d0*Qmean_n(i,j,5)+2.d0/3.d0*(Qmean(i,j,5)-dt*Qmean_n_k(i,j,5))     
    enddo  
    call comput_original_var
    call boundary
    
    
    
    
    

end subroutine 